WAVE COMPUTATION ON THE HYPERBOLIC DOUBLE DOUGHNUT 
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Abstract. We compute the waves propagating on the compact surface of constant negative curva- 
ture and genus 2. We adopt a variational approach using finite elements. We have to implement the 
action of the fuchsian group by suitable boundary conditions of periodic type. A spectral analysis of 
the wave allows to compute the spectrum and the eigenfunctions of the Laplace-Beltrami operator. 
We test the exponential decay due to a localized dumping and the ergodicity of the geodesic flow. 



The Hyperbolic Double Doughnut K is the compact manifold of negative constant curvature with 
two holes. We can define it by the quotient of the hyperbolic Poincare disc D, by some Fuchsian 
group r. Alternatively, we can construct it as the quotient of the so-called Dirichlet polygon, or 
fundamental domain ^ C D by a suitable relation of equivalence ~ : 



This beautiful object has many fascinating properties as regards the classical and quantum chaos 
(classical references are [1], [5]). Several important computational investigations of the spectrum 
were performed by using a stationnary method by R. Aurich and F. Steiner [4J. Moreover there 
has been much recent interest for the cosmological models with non trivial topology (a seminal 
work is the famous "Cosmic Topology" by M. Lachieze-Rey and J-P. Luminet [7:). In this context, 
K has been studied as a paradigm in [3], where a schema based on the finite differences on a 
euclidean grid was used to solve the DAlembertian. In this paper we compute the solutions of 
the wave equations in the time domain, by using a variational method and a discretization with 
finite elements on very fine meshes. The domain of calculus is the Dirichlet polygon, therefore 
the initial Cauchy problem on the manifold without booundary K, becomes a mixed problem on 
J- and the action of the Fuchsian group is expressed as boundary conditions on dT ^ analogous to 
periodic conditions. These boundary constraints are implemented in the choice of the basis of finite 
elements. By this way we obtain very accurate results on the transient waves. We test these results 
by performing a Fourier analysis of the transients waves that allows to find the first eigenvalues 
of the Laplace-Beltrami operator Ak on K. We compute also the solutions of the damped wave 
equation 



When < a £ L°°(K) and a > on dT^ the ergodicity of the geodesic flow assures that the 
geometric control condition of Ranch and Taylor [9 is satisfied. Our numerical experiments agree 
with their theoretical results, stating that the energy decays exponentially. 



In this part we describe the construction of the Hyperbolic Double Doughnut. First we recall 
some important properties of the 2-dimensional hyperbolic geometry. It is convenient to use the 
representation of the hyperbolic space by using the Poincare disc 



I. Introduction. 



K = D/r = Tj 



II. The Hyperbolic Double Doughnut. 




D := {{x,y)£R\ x' + y^ < l] 
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endowed with the metric expressed with the polar coordinates by 

("•2) = (1,^2)2 ^^' + 4(3^72)2 V = (i_Xy2)2 [d^' + dy']- 

It is useful to use the complex parametrization z = x + iy. We have to carefully distinguish the 
euclidean distance 

(II.3) d{z,z') =\ z- z' \, 



and the hyperbolic distance associated with the hyperbolic metric, given by 

2\z - z'? 

(11.4) COShdi/(z,2;') = 1 + (i_|^|2)(l_|^/|2) - 

We remark that 

d(0,z)=tanh^^^ 

hence the euclidean circles centerd in are hyperbolic circles, and more generally, all the hyperbolic 
circles {z'; dui^z' ^ zq) = i?}, with i? > 0, G D, are euclidean circles. The invariant measure dun 
on the Poincare disc allows to compute the area of any Lebesgue measurable subset X C D by the 
formula 

(11.5) m{^) = -|z|2)2 

in particular the hyperbolic area of a disc Dh{0, R) := {z; dniz, 0) < i?} is : 

(11.6) i^h{Dh{0,R)) = 47rsinh2 ( ^) . 



The group of the isometrics of D is generated by three kinds of transformations. 
(1) The Rotations of angle i^o S 1^ 

R^,iz) = e"^''z, 
and so R^^ is defined in the (x, y)-coordinates by the matrix 

R -(^'^^ M 



(2) The Boosts, or Mobius transforms, associated with tq G 



cosh + sinh ^ 
sinh ^z + cosh ^ 



2 ^ ' 2 
expressed in (x, y)-coordinates by the matrix 

rp _ ( cosh ^ sinh ^ 
""0 ~ \^ sinh^ coshf 

We remark that 
therefore 

(11.7) VZG]-1,1[, dH{z,Tr,{z)) = To. 

(3) The Symmetry 



cosh.dH{z.,TT-Q{z)) = 1 + 2j-_^ \„\2\2 
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Finally we recall that the geodesies of the Poincare disc are the diameters and all the arcs of 
circles that intersect orthogonally the boundary of the disc. 

Now we are ready to describe the double doughnut that is the quotient of the hyperbolic plane 
by the Fuchsian group of isometries, F, generated by the four transforms go^ 91: 92, 93, where 

(n.8) gk = Tr, R^ki 

with 



(II.9) tanh^ = yV2-l. 

The matrix of g^ is given by : 

These isometries gt satisfy the relation : 

(11.11) {9091 ^ 929^^) [go^ 9192^ g^) = Id- 
The Hyperbolic Double Doughnut is the quotient manifold 

(11.12) K:=D/F, 

endowed with the hyperbolic metric ds^ induced by ds^ . We know, see e.g. [1], [2], [5], that K is a 
two dimensional compact manifold, without boundary, its sectional curvature is constant, equal 
to —1, and its genus, that is the number of "holes", is 2. The geodesic flow is very chaotic : it is er- 
godic, mixing (theorems by G.Hedlung, E. Hopf), Anosov and Bernouillian (D. Ornstein, B. Weiss). 

To perform the computations of the waves on the doughnut, it is very useful to represent it by 
a minimal subset T <ZY) and a relation of equivalence ~ such that 

(11.13) K = J^/ ~ . 

When T is choosen as small as possible, it is called Fundamental Polygon. We take 

(11.14) .F:={zGD; Vi = 0,...,3, |5.(^)|>|z|, I^^H^)! > 1^1} • 

We can see that T \s & closed regular hyperbolic octogon, of which the boundary dJ- is the union 
of eight arcs of circle, that are parts of geodesies of D. We denote Pj, j G Zg, the tops of and 
PjPj+i the eight wedges. The action of F on the boundary is described by : 



(11.15) 



and 

P1P2 = 53(^6^5), P2P3 = 92{P7P6), 

^ ' P3P4 = giiPsPj), P4P5 = 90{PlP8)- 

We define the relation ~ by specifying the classes of equivalence z of any z G J^, z := F({2;}) n JF, 



Pi 


= 93{P6), 


P5 


= 93HP2) 


P2 


= 52(^7), 


Pe 


= 92\P^) 


P3 


= 9l{P8), 


P7 


= 9^\Pa) 


P4 


= 9o{Pi), 


Ps 


= 9o\P5) 



I.e. 



(11.17) z i = {z}, 

(11.18) Pj = {Pi, P2, Pa, P4, P5, P%, Pi, i^s}, 

(11.19) Z G dT\Pj ^ Z = {Z, Zequiv}, 
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Figure 1: The Fundamental Domain T . 



where according (11.16) 
(11.20) 



We give some metric relations. We denote P = PiPs H M^, P' = go{P) = Tj-^{P). We have 

dH{P,P')=2dH{0,P')=Tl, 

and so 



(11.21) 



PiPs c { 



z: \ X + 



'1 + V2 



^/2 - 1 



and by using the rotations we also have : 

dH{P^, Pi+l) =n= 2dH{P, Pi) = 2dH{0, P), 

dH{0,Pi) = T2 with tanhy = 2-3, d{0,Pi) = 2^^. 
Finally the area of the fundamental domain is ijlh{^) = 47r. 
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III. The Waves on the Doughnut. 
The Laplace Beltrami operator associated with a metric g, is defined by 

^ --d^g^^VWldu, g-' = {gn , IffI = I det <7^.|. 



\/\9\ 

We consider the lorentzian manifold x K endowed with the metric 

(111.1) g^ydx^dx" = dt^ - dsl^. 
We study the covariant wave equation 

(111.2) - AkV- = 0, 
and more generally, the damped wave equation 

(111.3) dfi; - Ak^ + a9tV = 0, 

whith < a G L°°(K). Here Ak is the Laplace Beltrami operator on K. Since K is a smooth 
compact manifold without boundary, K endowed with its natural domain {u € L^(K); Ax^t G 
L^(K)} is self- adjoint and the global Cauchy problem is well posed in the framework of the finite 
energy spaces. Given m G N, we introduce the Sobolev space 

(IIL4) H'^iK) := {u G L'^{K), VfjU £ L'^{K), \ a \< m} 

where Vj/ are the covariant derivatives. We can also interpret this space as the set of the distribu- 
tions u G HJ^^CD) such that uo g = u for any 5 G L. Then for all V'o £ -f^^(K), ipi G L^(K), there 
exists a unique V G C° {Rp, H^{K)) n (Mj+; ^^(K)) solution of ( |IIL3[ ) satisfying 

(IIL5) V(i = 0)=Vo, 5fV(t = 0) = Vi, 

and we have 

(IIL6) f I dti^{t) |2 + I Vui^it) dnH+ [ [ a \ dtip{t) p d^^Hdt = Cst. 

Jk Jo Jk 

To perform the numerical computation of this solution, we take the fundamental polygon as the 
domain of calculus. Then the Cauchy problem on x K is equivalent to the mixed problem 

(IIL7) dttiJ - ~ [ + dyyi>'\ + a{x, y)dtiP = 0, (t, x, y) G M+ x 

with the boundary conditions 

(IIL8) V(i, z) G M X z ^i;{t,z) = i;{t,z). 

We denote = {u G L^(J^), Va G N^, \a \< m, d'^yU G L'^{T)} the usual Sobolev space 

for the euclidean metric, and we introduce the spaces W^{J-) that correspond to the spaces 
//'"(K) : 

(IIL9) W'^iJ') := {u\r. u G i^;'^,(D), Vff G L, uog = u], 

endowed with the norm 

(in. 10) ||it||vF"' := ll'"|:'="ll_f/™(:F)- 

In particular, we have 

(III.ll) W^{T) = {ue H^{T), z ^ z ^ u{z) = u{z)] , 
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and for all V^p £ V F ^-^), V^i e there exists a unique ^ G C° {Rf; W^J')) D (Mi+; L'^{J')) 

solution of (III.5), ( |III.7 ) and ( |III.8 ), and we have the energy estimate : 



(III.12) 



Y ^ I dti^{t,x,y) p + I dxi'{t,x,y) ^ + | dy'4'{t,x,y) |^ dxdy 

[1 — x"^ — y^j^ 



+ / / 7^1 9^ — ^ I c^tV'Ct, X, y) p dxdy = Cst. 

Jo Jr (1 - a;^ - 



Since a G we have a result of regularity : when V'o £ ^^(•^)) V'l ^ then E 

C° (Mt;M^^(:r))nC^ (Mt;M^^(J^))nC^ (Mt;L2(jr)). In this case the mixed problem can be expressed 
as a variational problem : -0 is solution iff for all cf) G {J-) , we have : 

^W{t, z)(p{z)dxdy + — / ^ ^ip{t, z)(j){z)dxdy 



dt"^ (1 - x2 - y2)2 '^^ ' ^-^^ ^ « ' 7^ (1 - x2 - y2)2 

Ax,y'4'{t, z)<l){z)dxdy = 0. 



To invoke the Green formula, we denote v{z) the unit outgoing normal at z G dT. We suppose 
that PkPi = ga (PiPj). Then for z G P.P^- 

Since u o ga = u, we have for ti G W'^{T) 

du{z)U{z) = gaW{z)].Vu{ga{z)) = -d^(^g^(^))U{ga{z)) . 

We deduce that for u eW'^{T), v ^ W^{T), we have 

v{z)d^i^:,)ud\{z) = - v{z)d^(^^)udX{z), 
JPkPi 



and therefore 



v{z)di_,(^z)'^dX{z) = 0, / Ax,yu{z)v{z)dxdy = — dxudxV + dyudyvdxdy. 



We have proved the following 



Theorem III.l. Given i/jq G W'^{T), tpi G W^{J^), the solution if: of the Cauchy problem {III. 3), 
|///. 5p , is the unique solution satisfying (III. 5) of the variational prohl 



lem 



(III.13) 

d f 4a{z) 



WHJ^), ^ j^^^^:^j^i^{t,z)(t>{z)dxdy + 



jip{t, z)(i)[z)dxdy 



dt (1- I Z |2)2 

+ / dxtpit, z)dx(l}{z) + dyip{t, z)dy(j){z)dxdy = 0. 



We solve this variational problem by the usual way. We take a family V/j, < /i < /iq, of finite 
dimensional vector subspaces of W^{J-). We assume that 

(111.14) ^OKhKhoVh = W\J'). 

We choose sequences V'o,/i) ^ such 
We consider the solution ^ph G C'^iRt', Vh) of 
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(III. 15) 



[ 4 , / X , / X 7 , d [ Aaiz) , , X , / X 1 I 
^ J J, (1- I z |2)2 ^fc(^' zjcphizjdxdy + ^ (1- | z |2)2 ^^(^' zjcphizjdxdy 



+ 



^ dxi'hit, z)dx(t)h{z) + z)dy(t)h{z)dxdy = 0, 



satisfying ■0/j(O,.) = V'o,/i(-)) dttphi^,-) = V'i,fe(-)- Thanks to the conservation of the energy, this 
scheme is stable : 

Vr>0, sup sup + ll:nV'/iWllL2 < oo. 

o</i<?io o<t<r "I 



Moreover, when ijj e (M^; T^-'^(jr)), it is also converging : 



vr>o, sup \\Mt) - i^{t)\\w^ + - TM^^h' ^ ^^o. 

o<t<T at at 



If we take a basis ( e 



"'J^ ) oiVh, we expand tph on this basis : 



and we introduce 



Then the variational formulation is equivalent to 

(III. 16) MX" + BX' + KX = 0. 

This differential system is solved very simply by iteration by solving 

AT 

(III. 17) M(X"+^ - 2X" + X"-^) + — D(X"+^ - X"-^) + {ATfKX'' = 0. 

We know that this scheme is stable, and so convergent by the Lax theorem, when 
(III.18) 

and if there exists K > such that 
(III. 19) 



< KX, X > 4 
sup — — < 



the CFL condition 
(III.20) 



xJo<MX,X > Ar2' 
K 

V/t e]0,^o], V(?!)/i e Vh, \\Va:,y(l)h\\mj^) < 

KAT < 2h, 



1- 
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is sufficient to assure tlie stability and tlie convergence of our sclieme. 



IV. Numerical resolution 



IV. 1. Mesh. First of all we construct the boundary dJ- from the equation (11.21 ) and we perform 
a discretization that is equidistant for the hyperbolic metric. Next we use the mesh generators 
Emc2^^and bamg"'"'^created by INRIA. If we only use Emc2"'"'^, the mesh contains too many 
vertices and is not suitable for the hyperbolic metric. So a first mesh is created by Emc2'^'^. We 
also consider a circle which is uniformly discretized with the same hyperbolic step than the exterior 
geometry. The radius is choosen as the final mesh is almost uniform. At last, we impose on every 
point of the exterior and interior geometry a metric, in the sense of bamg"*"^. This software can 
next create a mesh which is more uniform, with respect to the hyperbolic metric, than the first 
mesh, and that has a reasonnable number of vertices. 

To test the uniformity of the mesh, we compute the extrema of the hyperbolic distance between 
two neighbor vertices. As a check of the accuracy of the meshes we evaluated the area of the 
polygons created by the meshes, and we compared to 47r (area of the domain). Here are some 
examples: 



Mesh number of vertices maxdn mindH area /An 

Meshl : 7448 0.087 0.027 1.00012 

Mesh2 : 17574 0.049 0.0177 1.00007 

Mesh3 : 37329 0.036 0.012 1.00003 

MesM : 67517 0.027 0.009 1.000018 



In our meshes, the greater hyperbolic distance between consecutive vertices is not reached near 
the exterior boundary. To give an idea of the accuracy of this dicretization, we show in the following 
figure, a very rough mesh : 



Wave Computation on the Hyperbolic Double Doughnut 



9 




Figure 2: A very rough mesh with 1756 vertices. 



IV. 2. Vh space. We construct the finite element spaces of Pi type. We note 7^ all triangles of a 
mesh, and J^h = ^KeT^K- 

Vh ■.= {v:Th^^,ve C\Th), yK € % v\K e ^i{K), M ^ M' ^ v{M) = v{M')] 



If Mi and Mj denote two vertices of the mesh, we define a basis [e -] oiVh by: 

V ^ Jl<j<Nh 

(1) If M, dJ' : e){M,) = 5ij 

(2) IfM,isaP,point:e,^(M.) = {; ^^=^ 

(3) If M.j e dJ^, and is not a R- point: e^fM^) = | I ^{,^* ^. 

^ ' ■' 1 ^ ' y [) otherwise 

In particular, we have to determine the equivalent points on dJ-. To that, we write a program 
implementing the relations (11.20). 

The number of nodes Nh is the sum of the number of the vertices which are not in dJ-^ the 
number of vertices which are on four consecutive arcs of dJ- without beeing a Pi point, and one 
(because all Pi points are equivalent to one of them) . 



IV. 3. Matrix form of the problem. and M(i, j) are found with a numerical integration 

using the value at the middle of the edges of the triangles. 
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The stiffness matrix IK and the mass matrix M are sparse and symetric matrices. So we choose 
a Morse stockage of their lower part, and all of the calculations will be performed with this stockage. 

To solve the linear problem we use a preconditionned conjugate gradient method. The precon- 
ditionner is an incomplete Choleski factorisation, and the starting point is the solution obtained 
with a diagonal preconditionner. 



IV.4. Initial data. We consider the case where the initial velocity = 0, i.e. = = 0, and 

we choose first initial datas with a more or less small support near a given point. For instance, for 
the wave depicted in Figure 3, we have taken 



(IV.l) tpoix,y) = lOOeioo^^+ioof^- 



0, for x^ + y'^> 



100' 



IV.5. Discretized energy. In order to see the stability of our method we perform Ed{t) the 
discretized energy at the time t: 



E{nAt) = 



At 



At 



It is well known that our schema is conservative when a = 0, hence must be invariant all along 
the resolution. We test this property with the previous initial data. 





EdiO) : 


EdilOO) : 








with At = 0.001 : 


with At = 0.0005 : 


Meshl : 


8455.89602005935 


8455.89602005865 




Mesh2 : 


8484.17400988788 


8484.17400988815 


8484.17400988936 


Mesh3 : 


8494.43409419468 


8494.43409419464 


8494.43409419511 


MesM : 


8498.39023175937 


8498.39023175945 


8498.39023175923 
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IV. 6. Eigenvalues. We test our scheme in tlie time domain, by looking for the eigenvalues of 
the hamiltonian when a = 0. Since the Laplace-Beltrami operator Ak on the Hyperbolic Double 
Doughnut is a non positive, self-adjoint elliptic operator on a compact manifold, its spectrum 
is a discrete set of eigenvalues —q^ < 0, and by the Hilbert-Schmidt theorem, there exists an 
orthonormal basis in L^(K), formed of eigenfunctions {^pg) C i7°°(K) associated to g^, i.e. 



(IV.2) 



(1 



y 



2\2 



We take Vo = Therefore any finite energy solution ip{t,x^y) of d^if) — AkV' = has an 

expansion of the form e*''*'0g(^) 2/) (such expansions exist also for the damped wave equation, 
when a > 0, see [6]). More precisely, if we denote <, > the scalar product in L'^(K), we write 

1 



^{t, X, y) = — {< dtij{0, .), 1 > t+ < i;{0, .),!>) 



(IV.3) 



An 



sin qt 
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To compute the eigenvalues q, we investigate the Fourier transform in time of the signal ip^t, x, y) 
in the case where dtip{0, x, y) = 0. We fix some large T >> 1, and we put ^^^{X) := 'ip{t, X)e^^^ dt. 
Then ~ CT'^{uP' — q'^)~^, T — > oo. Practically, during the time resolution of the equation 

we store the values of the solution at some points M, including the origin, P', Mq near P4, for the 
discrete time kAt, Ni < k < Nj. We choose the initial step Ni in order to the transient wave is sta- 
bilized, that to say NiAt is greater than the diameter of the doughnut, i.e. NiAt > 2 * dn^O, Pi) ~ 
4.8969. Then we compute a DFT of {iphikAt, M))Ni<k<Nf with the free FFT library fftw. Let us 



note (^'j(M))o, 



Nf-N,+ 



1 the result. If ^'^(M) := 



Ekio^'^'^N^ + k)At,x,y)e~^-f--'^^"'^\ 



we 



search the values jmaxi, jmax2, ■■■ for which (||^j(M)|| )j has a maximum. Then the eigenvalues 
found by the algoritm expressed as: 



2tt 



{Nf -Ni + l)At 



J max 



We have made a lot of tests by varying parameters such as: At, Ni, Nj, the mesh, the observation 
point M. With the initial data (IV.l ), we find the following values for q : 



1,96 ±0,02 2, 85 ±0,04 
4, 34 ±0,06 4, 83 ±0,06 
6, 00 ±0,06 6, 63 ±0,06 

These results agree with the results obtained with a stationnary method with a mesh of 3518 
vertices in 

Alternatively, we could also use the power spectrum and calculate the square of the modulus of 
'^Ux,y) 

1 



-IIL^(K) = ^ |< V'(0,.),1 >|^ 



sm 



2 (g+<^)T 



+ 



sm 



2 (c^-(?)r 



(w - qy 



+ 



(cos^ cjT — cos cjT cos qT) 



Therefore, for an eigenvalue q^: 



1^, 



'^lli2(K) ><V(0,.),Vg >' 



T 
2q 



IV. 7. Damped waves. We test our scheme for the damped wave equation (III. 3) when the damp- 



ing function a > is non zero (for deep theoretical results, see [S], [S], We know that the energy 
of any finite energy solution decays exponentially (uniformly with respect to the initial energy) iff 
the dumping a satisfies the assumption of geometric control introduced by J. Ranch and M. Taylor 
in [9]. This condition means 



(IV.4) 



a{x{t),y{t))dt = ±00 



for any geodesic (x(t),y(t)). Since the geodesic flow on the compact Riemaniann manifold with 
constant negative curvature is very chaotic (more precisely ergodic, mixing, Anosov, Bernouillian 
see e.g. [2], [5]), it is sufficient to have a > near dJ-. Nevertheless we constat an exponential 
decay for some solution, even if we choose a dumping function a equal to a positive constant on 



very small support that does not satisfy (IV.4| : a > only on one triangle and its close neighbors. 



The next figures are obtained with mesh3 and a defined by: 



a{x,y) = 0, for , I z |< 0.6 ; a(x,y) = 0.1, otherwise. 
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20 40 60 80 100 120 140 



Figure 4: Logarithm of the energy as a function of time. 




Figure 5: Solution at the origin. 
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